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Abstract 

Various Monte Carlo techniques are used to determine the complete phase diagrams of the square 
well model for the attractive ranges A = 1.15 and A = 1.25. The results for the latter case are in 
agreement with earlier Monte Carlo simulations for the fluid-fluid coexistence curve and yield new 
results for the liquidus-solidus lines. Our results for A = 1.15 are new. We find that the fluid-fluid 
critical point is metastable for both cases, with the case A = 1.25 being just below the threshold 
value for metastability. We compare our results with prior studies and with experimental results 
for the gamma-II crystallin. 
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I. INTRODUCTION 



It has been known since the the work of Gast et al^, and subsequently confirmed by 
others^, that the phase behavior of colloidal particles depends sensitively on the range of 
attraction between them. For sufficiently short range attractive interactions, the phase di- 
agram exhibits a solid-fluid coexistence curve that is subtended by a metastable fluid-fluid 
coexistence curve. Such phase behavior is also typical of certain globular proteins in solu- 
tion^. This has led scientists to model globular proteins in solution using the ideas developed 
in colloid science. Several models with short range attractive interactions 6 i 7 i 8 i 9 have been 
used to calculate the phase diagrams of various proteins using computer simulation. As 
well, many of these models have been shown to obey a kind of van der Waal's extended 
states behaviop^. These isotropic short-range models have also been used to characterize^ 
the nucleation rates for globular proteins. The simplest of these models is the square well 
potential, given by 

oo, r < a 

V ( r ) = { -e, o <r < Ad (1) 
0, r > Xa. 

The fluid- fluid coexistence curve for values of A > 1.25 is known 11 . However, to date, no 
direct simulation results of the complete phase diagram are available for A = 1.25 or for 
smaller ranges away from the adhesive-sphere limit. This choice of A = 1.25 for the the 
range of interaction is believed to be close to the threshold value below which the square 
well model becomes metastable. Indeed, theoretical and semi-analytical treatments confirm 
this hypothesis*^, although the threshold values for other short-range models is smaller. 
Also, the phase behavior at this value for the square well model has been compared to that 
of the gamma crystallin proteins^, mutants of which have been linked to genetic cataracts^. 
The square well model has also been used to determine nucleation rates at several ranges of 
A to better understand protein crystallization^ and is also of interest for comparing results 
with theoretical approaches modelling phase behavioi^, as well as for comparison with 
experimental studies. It is important, therefore, that these phase diagrams be accurately 
calculated, both at and below A = 1.25. In this paper, we present such results using standard 
Monte Carlo and parallel tempering techniques. 
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II. COMPUTATIONAL METHODS AND DETAILS 



A. Solid-fluid coexistence 

Kofke 17 showed that the solid-fluid coexistence curve can be obtained by solving the 
first-order Clausius-Clapeyron equation, given by^ 

dP _ Ah 

~df3 ~ ~JEv~ ^ ' 

where P is the pressure, /3 — 1/T is the inverse temperature, h is the molar enthalpy, and 
v is the molar volume, respectively. The method, known as the Gibbs-Duhem method, 
requires that two isobaric-isothermal (NPT) Monte Carlo simulations be carried out in 
parallel so that information can be used to calculate the next state and, thus, points along 
the whole coexistence curve can be obtained. This method has been employed for solid-fluid 
coexistence for a variety of models^^*- One caveat is that an initial coexistence point 
must first be known in order to integrate the Clausius-Clapeyron equation. To do this, the 
equation of state along an isotherm is obtained using NPT simulations. Once the equation 
of state is known, we fit the fluid-solid lines with an equation of the form 

PP = i ~ + Ky^Y + c{-^-f (3) 
1 — ap 1 — ap 1 — ap 

where p is the number density. The solid line is fit to a 2 nd order polynomial equation of 
the form ap 2 + bp + c. Integrating these two equations yields the chemical potentials for the 
liquid and solid, as shown in equations (4) and (5), respectively. 



pA 2 . bla — c/a 2 + 1 c/ 2a 2 + bp cp 2 

Pin = In / + -L— -L + -L -f + P 3 - (b a - c/2a 2 + 1 4 

' ap 1 — ap (1 — ap) z (1 — ap d ' 



(3p s = 2ap + b(\np + 1) - (ap* + b\np* - c/p*) + (3f ex (p*) + \nA 2 p* - 1 (5) 

It should be noted that eq.(5) requires the knowledge of the free-energy at the reference 
density p*. To calculate this, we use the Frenkel-Ladd^ method of coupling a solid to 
harmonic springs, referred to as an Einstein lattice. The benefit of coupling the solid in 
this way is that the free energy of the Einstein lattice can be easily calculated analytically, 
allowing one to obtain the free energy of the system of interest. We consider a system that 



is dependent upon a coupling parameter £ such that the total energy of the system may be 
written as 



N 



u(0 = u o + zu = u o + zJ2( 




(6) 



i=l 



where N is the number of particles, f 0i i the position of the lattice site to which the particle i 
is assigned, U a the energy of the system of interest, and r*j are the positions of the particles. 
As the coupling parameter £ becomes large, the system becomes more strongly coupled to 
the lattice. For very large values the system will behave as a non-interacting Einstein lattice. 
To verify this behavior, the mean-squared displacements of the system for various values of 
£ are computed and compared to those of a non-interacting Einstein lattice. Once this value 
imax is obtained, we can calculate the free energy of the square well model at a reference 
density p* by2i 



where the first term represents the free energy of the Einstein lattice. With the chemical 
potentials of both phases now known, it is straightforward to calculate a coexistence point 
and to subsequently perform the Gibbs-Duhem method, thereby obtaining the fluid-solid 
coexistence curve. 

To obtain a coexistence point, NPT simulations are first performed for N = 256 particles 
in a periodic simulation cell. In one NPT simulation, on average, one volume displacement 
is attempted for every N attempts at displacing a particle. This is done because a volume 
move is computationally more expensive than a particle displacement. Equilibration runs 
lasted for 10 million Monte Carlo steps and production runs ran for 20 million Monte Carlo 
steps, a step being an attempt at either displacing a particle or changing the volume of the 
cell, on average. In order to compare the chemical potentials for coexistence, we must first 
perform simulations on a system coupled to an Einstein lattice, as described previously. We 
perform regular Metropolis Monte Carlo simulations on a system of N = 256 particles with 
the constraint that the center of mass of the particles remains fixed^. If a particle is given 
a random displacement, all particles are subsequently shifted in the opposite direction to 
ensure that the center of mass is constant. We update the position of the center of mass 
every time a trial move is accepted. Thus, the shift in the center of mass is continually 





4 



m \ 

* r, 

500 1000 1500 2000 2500 3000 3500 4000 

FIG. 1: Comparison of mean squared displacements for a coupled solid with that of an Einstein 
lattice (solid line) for the range A = 1.25. 

updated in order to properly calculate the harmonic energy contribution 22 . Einstein lattice 
simulations were equilibrated for 5 million Monte Carlo steps and were run for a total of 
10 million Monte Carlo steps. These simulations are carried out at different values of the 
coupling parameter £. Fig. 1 shows our results for the mean squared displacements for 
A = 1.25. We use £ max = 4000 at A = 1.25 and ^ ma x = 3000 at A = 1.15 for the square 
well systems, respectively. The integral in eq. (5) is evaluated using a 10-point Gaussian 
quadrature^ 3 - to obtain the free energy of the square well system. Once this free energy is 
obtained, the chemical potentials are plotted versus pressure to obtain a coexistence point, 
given by the point at which the two curves intersect. 

The Gibbs-Duhem integration is performed using two NPT simulations in parallel. Sim- 
ulations were equilibrated for 20 million Monte Carlo steps and ran for 40 million Monte 
Carlo steps in production. Each NPT simulation was performed under the same conditions 
as described previously. The coexistence pressures were calculated using a simple predictor- 
corrector algorithm 23 . With the predicted/corrected pressure, NPT simulations are done 
in parallel to obtain the next coexistence point, and the process is repeated until the full 
fluid-solid coexistence curve is obtained. This procedure was implemented for both systems 
at A = 1.15 and A = 1.25. 
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FIG. 2: Equation of state of the square well model for A = 1.15. Also shown as the dashed lines 
are our fits according to equations in the text. Data was collected at the temperature T = 1.0. 









T 




















* 


* 








si 






























* 




ye' 


* 






























* 



0.2 0.4 0.6 0.8 1 1.2 

P 

FIG. 3: Equation of state of the square well model for A = 1.25. Also shown as the dashed lines 
are our fits according to equations in the text. Data was collected at the temperature T = 1.0. 
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FIG. 4: Determination of a coexistence point for A = 1.25. The point at which the curves cross 
satisfies the condition for coexistence, namely, equal chemical potentials at equal pressures. 

B. Fluid-fluid coexistence 

To calculate fluid- fluid coexistence of the square well model at A = 1.25, we use the 
Gibbs ensemble Monte Carlo method^. Two physically separated but thermodynamically 
connected simulation cells are used to calculate both the less dense and more dense fluid 
phases, respectively. This method circumvents the problem of an interface, which would 
hinder sampling, by altogether removing it. The two cells are allowed to exchange particles 
and both separately undergo volume displacements such that the total number of particles 
N = Ni + N2 and total volume V — V1 + V2 remain constant. Particles are also displaced 
within each cell according to the regular Metropolis method. Simulations for N = 512 
particles were conducted on a periodic simulation cell. Simulations were equilibrated for 50 
million steps and produced for 100 million steps. We required that the chemical potentials 
of both phases at each state point be equal to ensure coexistence had been reached. 

At shorter interaction ranges, A, it becomes increasingly difficult to sample phase space. 
The system becomes non-ergodic and standard techniques used to calculate fluid-fluid coex- 
istence are unable to obtain coexistence points. As A is decreased, the critical temperature 
T c is decreased, and the particles tend to 'stick' together. Stickiness*^ is a phenomenon 

7 




usually associated with the adhesive hard-sphere model 25 . It is not unusual, therefore, that 
for sufficiently short-range square-well potentials stickiness begins to manifest itself. (An 
adhesive sphere interaction can be derived from a square-well interaction in the limit where 
the well becomes infinitesimally narrow and infinitesimally deep.) 

To overcome this difficulty, we employ parallel-tempering 2 ^ Monte Carlo to speed up 
equilibration of the system. Parallel tempering allows local Monte Carlo simulations (repli- 
cas) to communicate and exchange information between each other. The benefit of this 
is that systems that were unable to properly sample phase space are able to do so. Par- 
allel tempering is quickly becoming a standard method of sampling systems that become 
trapped in local energy minim a 26 i 29 i 30 . We set up global Monte Carlo simulations using 
grand canonical Monte Carlo simulations in each replica, where each replica is allowed to 
exchange particle-configurations according to 



where %i is the state of the i th replica and \i is the reduced chemical potential fx/ksT. The 
potential used is actually //* = /i — 3 ln[A/a] (as in Ref^). In what follows, we denote /i* 
as /i for simplicity. To ensure detailed balance is obeyed, a replica is selected at random 
and tested to swap particle-configurations with a neighboring replica according to the above 
probability. Within each replica three trial moves are attempted: 1) insertion/deletion 
trial moves are attempted according to standard Monte Carlo as adapted for use in the 
grand canonical ensemble 28 , 2) particle displacements are employed according to the regular 
Metropolis method, and 3) cluster moves^i are attempted to break-up any small clusters 
that may have formed. The last trial move is similar to those of particle displacements, 
however, extra care must be taken so that detailed balance is observed. Thus, a cluster of 
particles is displaced such that no new clusters are formed and no old clusters are destroyed. 

Periodic boundary conditions are used for a simulation cell of size L = 8a. Simulations 
were equilibrated for 10 million steps and production runs lasted 100 million steps. At each 
state point we chose to run our replicas at a common temperature with each differing in 
chemical potential //. We used 6 to 7 replicas per global simulation, choosing the chemical 
potentials such that approximately 20% of particle exchanges were accepted. Once a density 
distribution was obtained in this way, coexistence points were determined using an equal- 




(8) 
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FIG. 5: Phase diagram of square well model for A = 1.15 using Monte Carlo methods described 
in the text. Shown are the fluid-solid (v) and metastable fluid-fluid (A) coexistence curves along 
with our fit (solid line) to the latter. 

area criterion according to 



0.5. 



(9) 



To facilitate our grand canonical simulations, we employ histogram-reweighting 32 after 
each set of simulations. This method allows information at one state point to be obtained 
from a neighboring one, such that 



exp[(// - v)pV - (g - f3)uV}P { L M (p,u) 
E exp[(/i' - y)pV - 03' - [3)uV)Pi M (p, it) 

III. RESULTS AND DISCUSSION 



(10) 



We calculated the phase diagrams of the square well model for both A = 1.15 and A = 1.25, 
respectively. Our equations of state for both interaction ranges, simulated at the isotherm 
T = 1.0, are shown in Figs. 2 and 3. Our fits to the data are also shown in these figures, 
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FIG. 6: Phase diagram of square well model for A = 1.25 using Monte Carlo methods described 
in the text. Shown are the fluid-solid (v) and metastable fluid-fluid (A) coexistence curves along 
with our fit (solid line) to the latter. The inset shows a close-up view of the phase diagram near 
the estimated critical point (■). 



TABLE I: Results from fit to fluid-fluid coexistence curves 



A 


T c 


Pc 


4>c 


A 


B 


1.15 


0.565(3) 


0.404(3) 


0.212 


0.55(2) 


0.89(1) 


1.25 


0.765(2) 


0.390(4) 


0.204 


0.50(2) 


0.82(1) 



as described previously in the text. We emphasize the importance of obtaining very good 
fits to the data, as one must have accurate values of the chemical potentials of the two 
phases, vis-a-vis eqs. (4) and (5). Errors in fitting the data, particularly for the liquid, lead 
to errors in the value for a coexistence point and, ultimately, in determining the fluid-solid 
coexistence phase boundaries. Care was taken such that deviations between the data and 
the fit were usually less than 1.0%, no data point deviated from the fit by more than 3.0% for 
either model. Certain data points were obtained by running longer simulations to minimize 
fluctuations. 

After obtaining good fits to our isotherms, we obtained coexistence points for the square 
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FIG. 7: Comparison of our Monte Carlo results (a) to the cell model results of Ref. ^ (v) for the 
range A = 1.25. Also plotted are our estimates of the critical temperature T/T c . The metastability 
gap is represented as a solid line. 



well models. This is shown in Fig. 4 where a coexistence point for A = 1.25 has been 
determined by comparing the chemical potential versus pressure. We then employed the 
Gibbs-Duhem method for each model to obtain the fluid-solid coexistence curve. Figs. 5 
and 6 show the phase diagrams for both interaction ranges studied. It can be seen that for 
A = 1.15 the metastable fluid-fluid coexistence curve is far below the liquidus line. This 
curve clearly indicates that the threshold for metastability lies above above this value for 
the interaction range. Other studie a 1 ^ 1 ? 4 have been done on the Yukawa and square well 
models for extremely short range interactions. There, a phenomenon of isostructural solid- 
solid coexistence is observed, with the appearance of a solid-solid critical point. We observe 
no such coexistence at A = 1.15, although we have not extended our simulations to high 
densities where such a phenomenon has been observed at lower ranges. Fig. 6 shows the 
phase diagram for A = 1.25, showing that the fluid- fluid coexistence curve is metastable. We 
note that at the liquidus line just above the critical temperature, T c , simulations exhibited 
large fluctuations in the density. We therefore ran longer simulations in this region than for 
other parts of the fluid-solid phase diagram to minimize fluctuations. 
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FIG. 8: Comparison of our Monte Carlo results for both A = 1.15 (□) and A = 1.25 (v), respectively, 
to the gamma-II crystallin (•). 

Also shown in Figs. 5 and 6 are our attempts to fit the metastable fluid-fluid coexistence 
curves to the equation p± — p c = A\T — T c \ ± B\T — T c \", where T c and p c are the critical 
temperature and density, respectively, and j3 = 0.3258 (Ref 35 ) is the Ising expononet. Table 
I shows our values for both A = 1.15 and A = 1.25. Our estimates compare well with other 
predictions^*^ of the critical point for these two ranges. The critical density p c is shifted 
toward a higher density for A = 1.15 as compared to that at A = 1.25, as shown in Table I, 
but is smaller than the critical density predicted in (RefSai). 

As noted in the introduction, many models capture the qualitative characteristics of 
protein phase diagrams. It is interesting, however, that the Yukawa, MLJ, and 2n-n Lennard 
Jones models all have been shown^ to become metastable at the interaction range R = 
A — 1 = 0.13 — 0.15. As noted, this is not, however, the case for the square well model. Our 
results confirm that R ~ 0.25, as previous*^ studies have predicted. 

In another study 9 of the square well model at A = 1.25, the authors obtain the phase 
diagram by a combination of Monte Carlo and extrapolation techniques^. We compare our 
phase diagram at A = 1.25 to theirs in Fig. 7, converting to the units T/T c and volume 
fraction 0/0 c (p/p c ). It is observed that both studies predict similar behavior for both the 
solidus line and metastable fluid-fluid coexistence curve. The liquidus lines, however, are 
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not compatible. We quantify this discrepancy by measuring the metastability gap, defined 

by 9 - 



T L -T C 



n s - n* s 
n* s — rii 



(11) 



where n» are phenomenological values used in their model. Using this definition we find the 
gap values 0.240 for A = 1.15 and .035 for A = 1.25 for the metastability gaps. The former 
predictions 9 - overestimate the metastability gap by approximately 50% and 60%, respectively. 
The fact that our results for the solidus line agree with the early study is probably due to the 
fact that their choice for the number of contacts a given particle has within its interaction 
range in their cell model approximation is consistent with the fee crystal structure for this 
model. 

Much interest has been focused on the gamma crystalling 36 , a group of monomeric eye 
lens proteins. Certain mutations^ of these proteins has been linked to genetic cataracts and 
protein crystallization^. These proteins have been well studied^ 3 - and their phase behavior 
are similar to those of the square well model, warranting a comparison of the two. Fig. 8 
shows our square well phase diagrams to those 9 of the gamma-II crystallins. Although we 
compare our data with the gamma-II crystallin, it has been showr4i that the data for the 
fluid-fluid coexistence curve of the entire gamma crystallin family lies on the same curve, 
within the scatter of the data. Our values for <p c , shown in Table I, are close to that 
of the value for the family of gamma crystallins. Both models qualitatively capture the 
characteristics of the gamma protein, with the model at the range A = 1.15 more closely 
resembling it. The case A = 1.15 yields better agreement with the liquidus and solidus lines 
than the model with A = 1.25. For both models, however, the fluid-fluid coexistence curve 
is not as broad as that of the gamma-II crystallin. Indeed, some studie o 38 ' 39 have attempted 
to remedy this by including anisotropy in the interactions among particles in the square 
well model, thereby broadening the fluid-fluid coexistence curves. One study 3 ^ averages out 
the rotational degrees of freedom to obtain an effective temperature-dependent isotropic 
potential that better approximates both the liquidus and metastable fluid-fluid coexistence 
curves of the gamma-IIIb crystallin. 
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IV. CONCLUSION 



The square well model for the interaction ranges A = 1.15 and A = 1.25 have been 
obtained using Monte Carlo simulations, including parallel tempering. We find that the 
latter interaction range is just below the threshold value for metastability of the fluid- 
fluid coexistence curve. Finite-size effects have not been taken into account for either the 
liquidus line or metastable fluid-fluid coexistence curve near the critical point where critical 
fluctuations could affect our results. However, the metastability gap is sufficiently large that 
we believe our conclusion about the metastable nature of the coexistence curve for A = 1.25 
is correct. We have also compared our Monte Carlo results for the complete phase diagram 
at A = 1.25 to another study^ and find that this study does not predict the correct liquidus 
line. The complete phase diagram at A = 1.15 was obtained by using standard Monte 
Carlo techniques, overcoming the problem of 'stickiness' by using parallel tempering for the 
fluid-fluid coexistence curve. A comparison of these two phase diagrams with the gamma-II 
crystallin has been made. Better qualitative agreement is obtained at A = 1.15 with the 
experimental results for the gamma-II crystallin, but there is no quantitative agreement. 

Isotropic fluids cannot be expected to quantitatively match experimental systems. Al- 
though certain characteristics of the phase diagrams of proteins can be captured by these 
simple models, it is unlikely that a correct solidus line can be obtained. This is due to 
anisotropic interactions and the (non spherical) ellipsoidal shape of certain globular pro- 
teins, such as the gamma crystalling. Models such as the square well model, or similar 
ones, alone, cannot account for effects such as anisotropy, solvent-solute interactions, or 
hydration effects commonly found in experimental systems and are, thus, too simple to be 
quantitatively accurate. However, the isotropic models do seem to provide a useful first 
approximation which can then be extended to include effects such as anisotropy and hy- 
drophobicity. 
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